<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">











<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <title>Math - The Commons Math User Guide - Numerical Analysis</title>
    <style type="text/css" media="all">
      @import url("../css/maven-base.css");
      @import url("../css/maven-theme.css");
      @import url("../css/site.css");
    </style>
    <link rel="stylesheet" href="../css/print.css" type="text/css" media="print" />
        <meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1" />
      </head>
  <body class="composite">
    <div id="banner">
                    <span id="bannerLeft">
    
            Commons Math User Guide
    
            </span>
                    <div class="clear">
        <hr/>
      </div>
    </div>
    <div id="breadcrumbs">
          
  

  
    
  
  
    
              <div class="xright">      
  

  
    
  
  
    
  </div>
      <div class="clear">
        <hr/>
      </div>
    </div>
    <div id="leftColumn">
      <div id="navcolumn">
           
  

  
    
  
  
    
                   <h5>User Guide</h5>
            <ul>
              
    <li class="none">
                    <a href="../userguide/index.html">Contents</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/overview.html">Overview</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/stat.html">Statistics</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/random.html">Data Generation</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/linear.html">Linear Algebra</a>
          </li>
              
    <li class="none">
              <strong>Numerical Analysis</strong>
        </li>
              
    <li class="none">
                    <a href="../userguide/special.html">Special Functions</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/utilities.html">Utilities</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/complex.html">Complex Numbers</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/distribution.html">Distributions</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/fraction.html">Fractions</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/transform.html">Transform Methods</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/geometry.html">3D Geometry</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/optimization.html">Optimization</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/ode.html">Ordinary Differential Equations</a>
          </li>
              
    <li class="none">
                    <a href="../userguide/genetics.html">Genetic Algorithms</a>
          </li>
          </ul>
                                           <a href="http://maven.apache.org/" title="Built by Maven" class="poweredBy">
            <img alt="Built by Maven" src="../images/logos/maven-feather.png"></img>
          </a>
                       
  

  
    
  
  
    
        </div>
    </div>
    <div id="bodyColumn">
      <div id="contentBox">
        <div class="section"><h2><a name="a4_Numerical_Analysis"></a>
4 Numerical Analysis</h2>
<div class="section"><h3><a name="a4.1_Overview"></a>
4.1 Overview</h3>
<p>
         The analysis package is the parent package for algorithms dealing with
         real-valued functions of one real variable. It contains dedicated sub-packages
         providing numerical root-finding, integration, and interpolation. It also
         contains a polynomials sub-package that considers polynomials with real
         coefficients as differentiable real functions.
        </p>
<p>
         Functions interfaces are intended to be implemented by user code to represent
         their domain problems. The algorithms provided by the library will then operate
         on these function to find their roots, or integrate them, or ... Functions can
         be multivariate or univariate, real vectorial or matrix valued, and they can be
         differentiable or not.
        </p>
<p>
          Possible future additions may include numerical differentiation.
        </p>
</div>
<div class="section"><h3><a name="a4.2_Root-finding"></a>
4.2 Root-finding</h3>
<p>
          A <a href="../apidocs/org/apache/commons/math/analysis/solvers/UnivariateRealSolver.html">
          org.apache.commons.math.analysis.solvers.UnivariateRealSolver.</a>

          provides the means to find roots of <a href="../apidocs/org/apache/commons/math/analysis/UnivariateRealFunction.html">univariate real-valued functions</a>
.
          A root is the value where the function takes the value 0.  Commons-Math
          includes implementations of the following root-finding algorithms: <ul><li><a href="../apidocs/org/apache/commons/math/analysis/solvers/BisectionSolver.html">
          Bisection</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/solvers/BrentSolver.html">
          Brent-Dekker</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/solvers/NewtonSolver.html">
          Newton's Method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/solvers/SecantSolver.html">
          Secant Method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/solvers/MullerSolver.html">
          Muller's Method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/solvers/LaguerreSolver.html">
          Laguerre's Method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/solvers/RidderSolver.html">
          Ridder's Method</a>
</li>
</ul>
</p>
<p>
          There are numerous non-obvious traps and pitfalls in root finding.
          First, the usual disclaimers due to the way real world computers
          calculate values apply.  If the computation of the function provides
          numerical instabilities, for example due to bit cancellation, the root
          finding algorithms may behave badly and fail to converge or even
          return bogus values. There will not necessarily be an indication that
          the computed root is way off the true value.  Secondly, the root finding
          problem itself may be inherently ill-conditioned.  There is a
           &quot;domain of indeterminacy&quot;, the interval for which the function has
          near zero absolute values around the true root,  which may be large.
          Even worse, small problems like roundoff error may cause the function
          value to &quot;numerically oscillate&quot; between negative and positive values.
          This may again result in roots way off the true value, without
          indication.  There is not much a generic algorithm can do if
          ill-conditioned problems are met.  A way around this is to transform
          the problem in order to get a better conditioned function.  Proper 
          selection of a root-finding algorithm and its configuration parameters
          requires knowledge of the analytical properties of the function under
          analysis and numerical analysis techniques.  Users are encouraged
          to consult a numerical analysis text (or a numerical analyst) when
          selecting and configuring a solver.
        </p>
<p>
          In order to use the root-finding features, first a solver object must
          be created.  It is encouraged that all solver object creation occurs
          via the <code>org.apache.commons.math.analysis.solvers.UnivariateRealSolverFactory</code>
          class.  <code>UnivariateRealSolverFactory</code> is a simple factory
          used to create all of the solver objects supported by Commons-Math.  
          The typical usage of <code>UnivariateRealSolverFactory</code>
          to create a solver object would be:</p>
<div class="source"><pre>UnivariateRealSolverFactory factory = UnivariateRealSolverFactory.newInstance();
UnivariateRealSolver solver = factory.newDefaultSolver();</pre>
</div>
<p>
          The solvers that can be instantiated via the 
          <code>UnivariateRealSolverFactory</code> are detailed below:
          <table class="bodyTable"><tr class="a"><th>Solver</th>
<th>Factory Method</th>
<th>Notes on Use</th>
</tr>
<tr class="b"><td>Bisection</td>
<td>newBisectionSolver</td>
<td><div>Root must be bracketted.</div><div>Linear, guaranteed convergence</div></td>
</tr>
<tr class="a"><td>Brent</td>
<td>newBrentSolver</td>
<td><div>Root must be bracketted.</div><div>Super-linear, guaranteed convergence</div></td>
</tr>
<tr class="b"><td>Newton</td>
<td>newNewtonSolver</td>
<td><div>Uses single value for initialization.</div><div>Super-linear, non-guaranteed convergence</div><div>Function must be differentiable</div></td>
</tr>
<tr class="a"><td>Secant</td>
<td>newSecantSolver</td>
<td><div>Root must be bracketted.</div><div>Super-linear, non-guaranteed convergence</div></td>
</tr>
<tr class="b"><td>Muller</td>
<td>newMullerSolver</td>
<td><div>Root must be bracketted.</div><div>We restrict ourselves to real valued functions, not complex ones</div></td>
</tr>
<tr class="a"><td>Laguerre</td>
<td>newLaguerreSolver</td>
<td><div>Root must be bracketted.</div><div>Function must be a polynomial</div></td>
</tr>
<tr class="b"><td>Ridder</td>
<td>newRidderSolver</td>
<td><div>Root must be bracketted.</div><div></div></td>
</tr>
</table>
</p>
<p>
          Using a solver object, roots of functions are easily found using the <code>solve</code>
          methods.  For a function <code>f</code>, and two domain values, <code>min</code> and
          <code>max</code>, <code>solve</code> computes a value <code>c</code> such that:
          <ul><li><code>f(c) = 0.0</code> (see &quot;function value accuracy&quot;)</li>
<li><code>min &lt;= c &lt;= max</code></li>
</ul>
</p>
<p>
          Typical usage:
        </p>
<div class="source"><pre>UnivariateRealFunction function = // some user defined function object
UnivariateRealSolverFactory factory = UnivariateRealSolverFactory.newInstance();
UnivariateRealSolver solver = factory.newBisectionSolver();
double c = solver.solve(function, 1.0, 5.0);</pre>
</div>
<p>
          The <code>BrentSolve</code> uses the Brent-Dekker algorithm which is
          fast and robust.  This algorithm is recommended for most users and  the 
          <code>BrentSolver</code> is the default solver provided by the 
          <code>UnivariateRealSolverFactory</code>.  If there are multiple roots
          in the interval, or there is a large domain of indeterminacy, the
          algorithm will converge to a random root in the interval without
          indication that there are problems.  Interestingly, the examined text
          book implementations all disagree in details of the convergence
          criteria.  Also each implementation had problems for one of the test
          cases, so the expressions had to be fudged further. Don't expect to
          get exactly the same root values as for other implementations of this
          algorithm.
        </p>
<p>
          The <code>SecantSolver</code> uses a variant of the well known secant
          algorithm.  It may be a bit faster than the Brent solver for a class
          of well-behaved functions.
        </p>
<p>
          The <code>BisectionSolver</code> is included for completeness and for
          establishing a fall back in cases of emergency.  The algorithm is
          simple, most likely bug free and guaranteed to converge even in very
          adverse circumstances which might cause other algorithms to
          malfunction.  The drawback is of course that it is also guaranteed
          to be slow.
        </p>
<p>
          The <code>UnivariateRealSolver</code> interface exposes many
          properties to control the convergence of a solver.  For the most part,
          these properties should not have to change from their default values
          to produce good results.  In the circumstances where changing these
          property values is needed, it is easily done through getter and setter
          methods on <code>UnivariateRealSolver</code>:
          <table class="bodyTable"><tr class="a"><th>Property</th>
<th>Methods</th>
<th>Purpose</th>
</tr>
<tr class="b"><td>Absolute accuracy</td>
<td><div>getAbsoluteAccuracy</div><div>resetAbsoluteAccuracy</div><div>setAbsoluteAccuracy</div></td>
<td>
                The Absolute Accuracy is (estimated) maximal difference between
                the computed root and the true root of the function.  This is
                what most people think of as &quot;accuracy&quot; intuitively.  The default
                value is chosen as a sane value for most real world problems,
                for roots in the range from -100 to +100.  For accurate
                computation of roots near zero, in the range form -0.0001 to
                 +0.0001, the value may be decreased.  For computing roots
                much larger in absolute value than 100, the default absolute
                accuracy may never be reached because the given relative
                accuracy is reached first.  
              </td>
</tr>
<tr class="a"><td>Relative accuracy</td>
<td><div>getRelativeAccuracy</div><div>resetRelativeAccuracy</div><div>setRelativeAccuracy</div></td>
<td>
                The Relative Accuracy is the maximal difference between the
                computed root and the true root, divided by the maximum of the
                absolute values of the numbers. This accuracy measurement is
                better suited for numerical calculations with computers, due to
                the way floating point numbers are represented.  The default
                value is chosen so that algorithms will get a result even for
                roots with large absolute values, even while it may be
                impossible to reach the given absolute accuracy.
              </td>
</tr>
<tr class="b"><td>Function value accuracy</td>
<td><div>getFunctionValueAccuracy</div><div>resetFunctionValueAccuracy</div><div>setFunctionValueAccuracy</div></td>
<td>
                This value is used by some algorithms in order to prevent
                numerical instabilities. If the function is evaluated to an
                absolute value smaller than the Function Value Accuracy, the
                algorithms assume they hit a root and return the value
                immediately.  The default value is a &quot;very small value&quot;.  If the
                goal is to get a near zero function value rather than an accurate
                root, computation may be sped up by setting this value
                appropriately.
              </td>
</tr>
<tr class="a"><td>Maximum iteration count</td>
<td><div>getMaximumIterationCount</div><div>resetMaximumIterationCount</div><div>setMaximumIterationCount</div></td>
<td>
                This is the maximal number of iterations the algorithm will try.
                If this number is exceeded, non-convergence is assumed and a
                <code>ConvergenceException</code> exception is thrown.  The
                default value is 100, which should be plenty, given that a
                bisection algorithm can't get any more accurate after 52 
                iterations because of the number of mantissa bits in a double
                precision floating point number. If a number of ill-conditioned
                problems is to be solved, this number can be decreased in order
                to avoid wasting time.
              </td>
</tr>
</table>
</p>
</div>
<div class="section"><h3><a name="a4.3_Interpolation"></a>
4.3 Interpolation</h3>
<p>
          A <a href="../apidocs/org/apache/commons/math/analysis/interpolation/UnivariateRealInterpolator.html">
          org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator</a>

          is used to find a univariate real-valued function <code>f</code> which
          for a given set of ordered pairs 
          (<code>x<sub>i</sub></code>,<code>y<sub>i</sub></code>) yields
          <code>f(x<sub>i</sub>)=y<sub>i</sub></code> to the best accuracy possible. The result
          is provided as an object implementing the <a href="../apidocs/org/apache/commons/math/analysis/UnivariateRealFunction.html">
          org.apache.commons.math.analysis.UnivariateRealFunction</a>
 interface. It can therefore
          be evaluated at any point, including point not belonging to the original set.
          Currently, only an interpolator for generating natural cubic splines and a polynomial
          interpolator are available.  There is no interpolator factory, mainly because the
          interpolation algorithm is more determined by the kind of the interpolated function
          rather than the set of points to interpolate.
          There aren't currently any accuracy controls either, as interpolation
          accuracy is in general determined by the algorithm. 
        </p>
<p>Typical usage:</p>
<div class="source"><pre>double x[] = { 0.0, 1.0, 2.0 };
double y[] = { 1.0, -1.0, 2.0);
UnivariateRealInterpolator interpolator = new SplineInterpolator();
UnivariateRealFunction function = interpolator.interpolate(x, y);
double interpolationX = 0.5;
double interpolatedY = function.evaluate(x);
System.out println(&quot;f(&quot; + interpolationX + &quot;) = &quot; + interpolatedY);</pre>
</div>
<p>
          A natural cubic spline is a function consisting of a polynomial of
          third degree for each subinterval determined by the x-coordinates of the
          interpolated points.  A function interpolating <code>N</code>
          value pairs consists of <code>N-1</code> polynomials. The function
          is continuous, smooth and can be differentiated twice.  The second
          derivative is continuous but not smooth.  The x values passed to the
          interpolator must be ordered in ascending order.  It is not valid to
          evaluate the function for values outside the range 
          <code>x<sub>0</sub></code>..<code>x<sub>N</sub></code>.
        </p>
<p>
          The polynomial function returned by the Neville's algorithm is a single
          polynomial guaranteed to pass exactly through the interpolation points.
          The degree of the polynomial is the number of points minus 1 (i.e. the
          interpolation polynomial for a three points set will be a quadratic
          polynomial). Despite the fact the interpolating polynomials is a perfect
          approximation of a function at interpolation points, it may be a loose
          approximation between the points. Due to <a href="http://en.wikipedia.org/wiki/Runge's_phenomenon" class="externalLink">Runge's phenomenom</a>

          the error can get worse as the degree of the polynomial increases, so
          adding more points does not always lead to a better interpolation.
        </p>
<p>
          Loess (or Lowess) interpolation is a robust interpolation useful for
          smoothing univariate scaterplots. It has been described by William
          Cleveland in his 1979 seminal paper <a href="http://www.math.tau.ac.il/~yekutiel/MA%20seminar/Cleveland%201979.pdf" class="externalLink">Robust
          Locally Weighted Regression and Smoothing Scatterplots</a>
. This kind of
          interpolation is computationally intensive but robust.
        </p>
</div>
<div class="section"><h3><a name="a4.4_Integration"></a>
4.4 Integration</h3>
<p>
          A <a href="../apidocs/org/apache/commons/math/analysis/integration/UnivariateRealIntegrator.html">
          org.apache.commons.math.analysis.integration.UnivariateRealIntegrator.</a>

          provides the means to numerically integrate <a href="../apidocs/org/apache/commons/math/analysis/UnivariateRealFunction.html">univariate real-valued functions</a>
.
          Commons-Math includes implementations of the following integration algorithms: <ul><li><a href="../apidocs/org/apache/commons/math/analysis/integration/RombergIntegrator.html">
          Romberg's method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/integration/SimpsonIntegrator.html">
          Simpson's method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.html">
          trapezoid method</a>
</li>
<li><a href="../apidocs/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.html">
          Legendre-Gauss method</a>
</li>
</ul>
</p>
</div>
<div class="section"><h3><a name="a4.5_Polynomials"></a>
4.5 Polynomials</h3>
<p>
          The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/package-summary.html">
          org.apache.commons.math.analysis.polynomials</a>
 package provides real coefficients
          polynomials.
        </p>
<p>
          The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/PolynomialFunction.html">
          org.apache.commons.math.analysis.polynomials.PolynomialFunction</a>
 class is the most general
          one, using traditional coefficients arrays. The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.html">
          org.apache.commons.math.analysis.polynomials.PolynomialsUtils</a>
 utility class provides static
          factory methods to build Chebyshev, Hermite, Lagrange and Legendre polynomials. Coefficients
          are computed using exact fractions so these factory methods can build polynomials up to any degree.
        </p>
</div>
</div>

      </div>
    </div>
    <div class="clear">
      <hr/>
    </div>
    <div id="footer">
      <div class="xright">&#169;  
          2003-2009
    
          
  

  
    
  
  
    
  </div>
      <div class="clear">
        <hr/>
      </div>
    </div>
  </body>
</html>
